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£3 ; Abstract 
< 

We describe the computation of the Bondi news for gravitational radiation. 
We have implemented a computer code for this problem. We discuss the 
\ theory behind it as well as the results of validation tests. Our approach 

| uses the compactified null cone formalism, with the computational domain 

extending to future null infinity and with a worldtube as inner boundary. We 
calculate the appropriate full Einstein equations in computational eth form in 
Q^ . (a) the interior of the computational domain and (b) on the inner boundary. 

At future null infinity, we transform the computed data into standard Bondi 
coordinates and so are able to express the news in terms of its standard iV_|_ 
^Jy and iV x polarization components. The resulting code is stable and second- 

order convergent. It runs successfully even in the highly nonlinear case, and 
has been tested with the news as high as 400, which represents a gravitational 
radiation power of about W 13 M Q /sec. 
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I. INTRODUCTION 



In a series of papers in the 1960s Bondi and others fl|-|5| developed the theory of grav- 
itational radiation at null infinity, and this formalism has now become accepted practice: 
theoretical results about gravitational radiation at future null infinity (X + ) are usually ex- 
pressed in terms of the Bondi news function N. However, until now it has not been possible 
to calculate the news numerically, except in spacetimes with special symmetries. This paper 
describes the theory, as well as a code, for computing the news in general (asymptotically 
flat) spacetimes and in regimes that include the highly nonlinear case. The code runs with 
the news N as large as 400 (in geometric units in which N is dimensionless) . This is enor- 
mous and means that the code can cope with a power output of order N 2 c 5 /G ~ 10 60 Vt^ 
in conventional units. This exceeds the power that would be produced if, in 1 second, the 
whole Galaxy were to become gravitational radiation. 
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Most work in numerical relativity is done in the Cauchy "3 + 1" formalism, with the 
gravitational radiation estimated by perturbative methods [§-§|; these methods have not 
been tested on high-powered waveforms. However, here we are using the characteristic, 
or null cone, formalism with radial distance compactified so that X + is contained in the 
(finite) grid. The theoretical foundations for this approach were laid in the 1980s |9|-|P2" 



Numerically, it has been implemented for the model problem of the scalar wave equation fl4 
and in general relativity under various restrictive conditions [I5-IB|. 

In a previous paper |18j we were able to compute the news in the quasi-spherical case, 
i.e. with the nonlinear terms in Einstein's equations ignored. The notation and approach 
here follow that in jTj|. The computational domain is bounded by an inner (timelike or 
null) worldtube, and the outer boundary of the domain is X + . The computational domain is 
foliated into a sequence of null cones with initial data given on the first cone. Appropriate 
data on the inner boundary then enable us to use Einstein's equations to find the spacetime 
geometry throughout the computational domain (see figure 1 in |TB[). In addition to its 
other contributions, this work is also an important step towards the full implementation of 
Cauchy-characteristic matching as envisaged in ]Tj| (see also [p!9|-p3|), since the characteristic 
code is now complete. 

The first step in this paper is the calculation of the nonlinear terms in the Einstein 
equations for the Bondi-Sachs metric, and the expression of these terms in computational 
eth form [ 24, 25f] . This is done computationally using a symbolic programming language, and 
the results are described in Section J|. The inner boundary of the computational domain 
may be a timelike, or an incoming null, worldtube, and Section |T| discusses the required 
boundary data 0. Next, in Section [TV], we investigate behavior at the outer boundary 
X + . We use null cone coordinates that are defined quite generally in terms of a 2 + 1 
decomposition of the inner world tube. Consequently, they do not form an inertial Bondi 
system at X + . However, by carrying out a transformation to an inertial frame, we express 
the news in terms of the two standard polarization modes iV + and N x [2B|. Using computer 
algebra, these expressions are translated into computational eth form. 

At this stage of the paper, the theoretical background for the problem has been com- 
pletely specified. Section [V] considers various details of the numerical implementation, in 
particular those aspects that did not occur in the quasispherical case. We go on in Section |V| 
to describe the results of testing the code for stability and accuracy and demonstrate that 
the code is second-order convergent. 

We also perform a series of runs for an incoming asymmetric gravitational pulse incident 
on a Schwarzschild black hole. This is the fully nonlinear version of the classic black hole 
scattering problem first worked out in the perturbative regime |Z7| . As the magnitude of the 
pulse is ramped up various nonlinear effects become evident. The largest pulse investigated 
leads to a peak news of N = 400 for the gravitational radiation backscattered to X + . 

Some of the computational eth expressions found in this paper are quite long, and they 
are given in Appendices: Appendix A for the nonlinear terms in the Einstein equations, and 
Appendix B for the news. 
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II. THE NULL CONE FORMALISM 



First we summarize our notation and previous results about the Einstein equations in 
the null cone formalism. These earlier results were complete only in the quasispherical 
approximation |0| . We next proceed to discuss the computer algebra techniques that were 
used to extend these results to the full nonlinear case in computational eth form ||24| , and 
we present the completed results. 

A. Previous results 

We use coordinates based upon a family of outgoing null hypersurfaces. We let u label 
these hypersurfaces, x A (A = 2,3), label the null rays and r be a surface area coordinate. 
In the resulting x a = (u, r, x A ) coordinates, the metric takes the Bondi-Sachs form 

ds 2 = - ^e 2/3 (l + y) - r 2 h AB U A U B ^j du 2 - 2e 2p dudr - 2r 2 h AB U B dudx A 

+ r 2 h A Bdx A dx B , (1) 

where W is related to the more usual Bondi-Sachs variable V by V = r + W; and where 
h AB hsc — $c an d det{h AB ) = det(q AB ), with q AB a unit sphere metric. In analyzing the 
Einstein equations, we also use the intermediate variable 

Qa = r 2 e- 2 ?h AB U B . (2) 

We work in stereographic coordinates x A = (q, p) for which the unit sphere metric is 

q AB dx A dx B = -p^(dq 2 + dp 2 ), (3) 

where 

P=l + q 2 +p 2 . (4) 

We express q AB in terms of a complex dyad q A (satisfying q A q A = 0, q A q A = 2, q A = q AB q B , 
with q AB q B c = Sq and q AB = | (q A q B + (}aQb))- We fix the dyad by the explicit choice 

9 A = f(l,<) (5) 



with % = \J — 1. Note that we have departed from other conventions p5f to avoid factors of 



y/2 which are awkward in numerical work. For an arbitrary Bondi-Sachs metric, h AB can 
then be represented by its dyad component J = h AB q A q B /2, with the spherically symmetric 
case characterized by J = 0. The full nonlinear h AB is uniquely determined by J, since the 
determinant condition implies that the remaining dyad component K = h AB q A q B /2 satisfies 
1 = K 2 — J J. We also introduce spin- weighted fields U = U A q A and Q = QaQ A , as well 
as the (complex differential) eth operators 9 and 8. Tensor quantities and spin-weighted 
quantities are related as follows: 
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h 22 = —(J + J + 2K), h 23 



-- h 32 
P 



-^(J-J)i, h 33 = ^(2K-J-J), 



U 2 = j(U + U), = j{U - U)i, Q 2 = £±9., Q 3 



p 



(Q-Q)- 



(6) 



Refer to [p4] , p!q| for further details. 

The Einstein equations G^ v = decompose into hypersurface equations, evolution equa- 
tions and conservation laws. In writing the field equations, we follow the formalism given 
in 



TT1,|T2|| . We find for the hypersurface equations: 



A. 



W, 



Np, 

r -2 e 2?Q + Nui 

-r 2 (SJ + QK) tf + 2r 4 9 (r" 2 /5) p + N Q , 

K 2f3 TZ - 1 - e^M + K- 2 (r 4 (QU + W)) r + N w , 



where [23| 



n = 2K- mx + -m 2 j + s 2 J) + -^(Bjsj - Sjsj). 

2 V ' AK y ' 



(7) 
(8) 
(9) 

(10) 



(11) 



Next, the evolution equation takes the form 
2M, ur -(r-V(rJ) ir ) 



1 (r 2 W) r + 2r- 1 e p d 2 e fS - (r -1 W) r J + Nj. 



(12) 



The remaining independent equations are the conservation conditions, which are discussed 
in Section |T|. 

In the above equations Np, Njj, Nq, Ny/ and Nj represent the nonlinear aspherical terms 



(in a sense specified in QT8 1 ) . Expressions for these terms are known in the form of 
tensors and covariant derivatives, rather than complex spin-weighted fields and 9. We now 
calculate those nonlinear terms in the desired form. 



B. New results 



We have used computer algebra to calculate the nonlinear terms in two separate ways. 
In the first approach we start from the expressions reported previously fI5 |. In the second 
approach we start from the text-book definition of the Ricci tensor 



R 



fi pa Q pet I -pa -per pet per 

U Oi L [ii/ u u L an ~r 1 rra L iw 1 rru 1 n 



era fiu 



av iia' 



(13) 



We have checked that the two approaches gave the same results by substituting the total 
expressions (including the quasispherical and nonlinear parts) obtained for /3 r , U >r , Q >r , W )T 
and J iUr into the relevant part of the Ricci tensor, and confirming that the difference between 
the two results was zero. 
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The translation of a formula from tensor to eth formalism is ideal for computer algebra: 
it is straightforward and algorithmic, but lengthy. The first step is to expand any covariant 
derivatives in terms of partial derivatives and connection terms. Then using equation (^j), 
each h^B is expressed in terms of J, J and K, and the U A and Qa are expressed in terms of U 
and U, and Q and Q, respectively. At this stage everything is in terms of spin-weighted fields, 
but angular derivatives of these fields are partial derivatives with respect to the coordinates. 
The next step is to transform these partial derivatives to S and 8 operators. We use simple 
linear algebra applied to earlier results p4[ to obtain 



QA + QA-2ipsA . iffiA - QA + 2qsA) 
A , q = p , A iP = ^ '- (14) 

where A is a spin-weighted field with spin-weight s. For convenient reference we list all the 
spin-weighted fields that are used in the Bondi-Sachs metric and their spin-weights: 

Spin- weight: -2 -1 12 

Field: J Q, U P,K,W Q,U J 

Note further that the operator 9 increases the spin- weight by 1 and the operator 9 decreases 
the spin-weight by 1. For example, 9/3, 9J and d 2 U all have spin-weight 1. 

For quantities that are not spin-weight 0, the expressions for A <q and A tP involve q and p. 
However, the eth representation of geometric spin-weighted quantities (e.g. R r Aq A ) should 
not involve q and p, so if they appear due to one partial derivative they should be cancelled 
out by something else in the expression; also P should not appear in the final representation. 
These features provide useful checks in computer algebra calculations. 

The procedures described above are complete for quantities that involve one angular 
derivative, including second derivatives involving one angular derivative and one derivative 
with respect to u or r. Quantities involving two angular derivatives are unpacked one 
derivative at a time, as illustrated in the example 

_ , w , ( m + m ^ p((aw% + (5wq, p ) - 2 P (dw + mv) 

_ i(mW - d 2 W + 2qdW + B 2 W - mW - 2qdW) - 2p(cW + cW) ^ 

We present the results of our computer algebra calculations in Appendix 0. Then the 
complete set of null cone hypersurface and evolution equations is given by equations (0) 
to fljjD , augmented by equations ( |A1~1 ) to (|A6|) . These equations are suitable for numerical 



implementation in computational eth form. 

Most of the computer algebra was done using Maple, with some developmental work 
and checking on Mathematica and REDUCE. The Maple scripts are available on the World 
Wide Web (|http: / / artemis.phyast.pitt.edu/hpgnewsj) 



III. CONSERVATION CONDITIONS 

Our discussion up to this point has only addressed the six hypersurface and evolution 
equations, designated by Bondi as the "main" Einstein equations |IJ. They correspond to 
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the six components of the Ricci tensor itemized in Appendix [A]. Together with the equations 
R r a = 0, they form a complete set of components of the vacuum Einstein equations. Given 
that the main equations are satisfied, Bondi used the Bianchi identities to show that one of 
these, the "trivial" equation R r r = 0, is automatically satisfied. He further showed that the 
remaining three equations 

K = 0, q A R r A = 0. (16) 

are satisfied on a complete outgoing null cone if they are satisfied on a single spherical 
cross-section. By choosing this sphere to be at infinity, he identified these three equations 
as conservation conditions for energy and angular momentum. 

In the context of our present work, we obtain a solution by requiring that these conserva- 
tion conditions be satisfied on a worldtube T forming the inner boundary of the characteristic 
domain U . These conditions are automatically satisfied when matching across T to an inte- 
rior solution generated, say, by Cauchy evolution. Matching is anticipated to be a primary 
application of the characteristic algorithm. However, in carrying out a purely characteristic 
evolution of the exterior, based upon data on T and on an initial null hypersurface, these 
conservation conditions must be imposed separately. This applies to the stand-alone tests of 
the characteristic code presented in Section [VT|. An alternative is to take the limit in which 
the worldtube collapses to a world line, in which case the conservation conditions reduce to 
regularity conditions at the vertices of the null cones. This approach was used in numeri- 
cal evolution with an axisymmetric characteristic code [[H] but would be computationally 
prohibitive in the 3-dimensional case. 

The conservation conditions (|1(J) dictate which metric variables can be given freely on T 
and which are subject to constraints. They contain no second r-derivatives of the metric and 
are the analogues, with respect to an r-foliation, of the traditional momentum constraints 
in a Cauchy formalism. We can therefore express them as 

D i (K ij ->y ij K) = Q (17) 

in terms of the intrinsic metric 7^, the associated covariant derivative Di and the extrinsic 
curvature in the induced coordinates x % = (it, x A ) of the worldtube Y. 
From (P, The intrinsic metric is 

- fij dx i dx j = -e w ^-du 2 + r 2 h AB (dx A - U A du)(dx B - U B du). (18) 

In analogue to the 3 + 1 decomposition of the Cauchy formalism, a 2 + 1 decomposition of 
the timelike worldtube geometry leads to the identification of qab — r<2 hAB as the metric of 
the 2-surfaces of constant it which foliate the world-tube, e 2/3 V/r as the square of the lapse 
function and (— U A ) as the shift vector. The extrinsic curvature of Y is 

K iS = 6?5$V a n p (19) 

in terms of the unit normal n a = (r /V)^e^W a r. 

Part of the worldtube data is specification of the gauge. One simple choice is to set the 
shift to zero {U A = 0), and the lapse to one (V/r = e _2/3 ). In this case, the components of 
the extrinsic curvature are 
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K uu — 2/5 u 4~ ^ Quu,n (20) 



K uA = P,a + ^e- 2(3 g uAtr , (21) 
Kab = \{-9ab,u + e- 2 Pg A B,r), (22) 



with 



K = -2/?„ + -e~ 213 - \e-^g uu ^ (23) 
r 2 

where g uu , g u A and are obtained from ([[]). It is then easy to work out the specific form 
of the conservation conditions using the identity 

DiKy = JLzdjiy/^Kj) + \K mnl m \ t , (24) 

which holds for any symmetric tensor. 

In order to elucidate the content of these equations we lump under the symbol K.p or K A 
any purely null-hyper surf ace term (composed out of (3, U A , V and Kab and their r and x A 
derivatives) together with any purely worldtube term (composed out of Uab and its u and 
x A derivatives). Then the ti-component of ( K7h takes the form 

A« = ^; ( 25 ) 

and the x A -components are 

(e- 2/3 ^A,,-2/?, A ), M = /C A , (26) 
which may be re-expressed in computational eth form as 

Q,u = -2$P,u - Q A ICa. (27) 

Equations ( p5|) and fl2"TD are generalizations to a finite worldtube of the energy and 
angular momentum conservation laws found by Bondi. They imply that (3 and Q cannot be 
given freely on the the worldtube but are subject to constraints. Since the worldtube values 
of j3 and Q enter as integration constants in the hypersurface equations, the evolution scheme 
must either update their values in accord with equations ( |2"oD and (E7|) or prescribe them 
by matching to an interior solution. (In the test cases presented in Section [VI] we match 
to an analytic interior solution). The remaining integration constants can then be posed 
freely on the worldtube. This gives the formal result (in our present choice of gauge) that 
specification of Iiab on the initial null hypersurface and on the worldtube and specification 
of 13 and Q on the initial slice of the worldtube determine the future evolution of the system 
in the region outside the worldtube. 

A similar nullcone-worldtube evolution scheme holds for a general choice of lapse and 
shift, although the detailed equations are considerably more complicated. An example where 
the worldtube is null is presented in Section |V11| . 
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IV. NULL INFINITY 



The spacetime can be conformally compactified in terms of the metric ds 2 = £ 2 ds 2 , where 
I — 1/r. In (u,£, x A ) coordinates, it takes the form 

ds 2 = - (e 2/ W - h AB U A U B ) du 2 + 2e 2l3 dud£ - 2h AB U B dudx A + h AB dx A dx B . (28) 

Here £ is a conformal factor with future null infinity X + given by £ = 0. The physical 
quantities governing the total energy and radiation power from the system are constructed 
from the leading coefficients (functions of u and x A ) in an expansion of the metric in powers 
of £, in accord with asymptotic flatness at X + . The coefficients H, H AB , cab and L A which 
enter the construction of the news function are defined by the expansions, (3 — H + 0(£ 2 ), 
h AB = H AB + £c AB + 0(£ 2 ), U A = L A + 2£e 2H H AB D B H + 0(£ 2 ) and £ 2 V = D A L A + 
£e 2H (TZ/2 + 2D A D A H + AD A HD A H) + 0{£ 2 ). (Note here that we have related some of the 
expansion coefficients by using the asymptotic content of the Einstein equations). 



A. The News 

The gauge freedom in the conformal metric ds 2 is fixed by the gauge conditions adopted 
on the inner worldtube V. In order to construct the news function, it is convenient to refer 
to a conformal Bondi frame || with metric ds 2 = Q 2 ds 2 = uj 2 ds 2 , with Q = uj£, which 
satisfies the gauge requirements that Q AB := g AB \x+ = u 2 H AB is intrinsically a unit sphere 
metric at X + and that ^ Q/3 V a f2V/3fi = 0(Q 2 ). The conformal Bondi frame corresponds to 



an asymptotically Minkowskian inertial frame. Reference [131 discuses the calculation of the 
news in an arbitrary conformal frame relative to that in a Bondi frame. 

X + is geometrically a null hypersurface with null vector tangent to its generators given by 
n a = ^ a ^V p£l\x+ or ; equivalently, n a = g al3 V p£\x+ > with n a = uj^h ". In order to complete 
a basis for tangent vectors to X + , let Q a be a complex field in the neighborhood of X + 
satisfying Q a V ' a Vt = 0(0,), ~g a pQ a Q fi = 0(0) and g aP Q a Q p = 2 + 0(0). In the conformal 
Bondi frame, the news function can then be expressed in the form |28| 

N = ^Q a Q p Vj7 p O (29) 

evaluated in the limit of X + . (Our conventions are chosen so that the news reduces to Bondi's 
original expression in the axisymmetric case 0]). In terms of the g a p frame, with conformal 
factor Q = Q/u = £, we then have 

N = l -Q a Q\^^- - cV.V,-). (30) 

We can determine uj on X + in the g a frame by solving the elliptic equation governing 
the conformal transformation of the curvature scalar intrinsic to the 2-geometry of a u = 
constant slice of X + , 

1I = 2(uj 2 + H AB D A D B \oguj), (31) 
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where 1Z is defined in equation ([11]). The condition that g^V a QV pQ = 0(Q 2 ) is consistent 
with setting dgu = but gives a relation for the time dependence of u. Noting that 
g ap V a iVpi = e- 2H D A L A £ + 0(£ 2 ) and that g aP V pi = n a + 0(l), we obtain 2n a W a hg lu = 
— e~ 2H D A L A . This is used to evolve uj given a solution of (|3l|) as initial condition. 

In order to obtain an explicit expression for the news (|30| ) in the g a p frame we need a 
representation of Q*. The choice is independent of the freedom Q 13 — » Q 13 + / yn 13 , which 
leaves ( j5H ) invariant. However, it is important for physical interpretation to choose the 
spin rotation freedom Q 13 — > e~ %OL Q^ to satisfy n a 'V a Q 13 = 0(Q), so that the polarization 
frame is parallely propagated along the generators of X + . This fixes the polarization modes 
determined by the real and imaginary parts of the news to correspond to those of inertial 
observers at X + . 

We accomplish this by introducing the dyad decomposition h AB = {m A fn B + fn A m B )/2 
and fixing the spin rotation freedom of m A by setting 

1/4 



J\ 'La. KK+1 




mA =[j) ( V ' \/ " <f ' ' I / ^77=^7 ) ■ ( 32 ) 

so that H AB = (M A M B + M A M B )/2 where M A = m A + 0(Q). Here, in order to make 
m A independent of the phase of q A (the basis vector used in the 9 computations), we have 
introduced the phase factor e* 7 = (J/J) 1//4 . Now, setting Q 13 = e~ m u;~ 1 (0, 0, M A ) +777,^, 
the requirement of an inertial polarization frame, n a 'V a Q 13 = 0(Q), determines the time 
dependence of the phase a. We obtain 

2{d u + L A d A ){ia + logcu) = +H AC M c {d u M A + L B d B M A - M B d B L A ) (33) 

or, eliminating the time derivative of to, 

2i(d u + L A d A )a = D A L A + H AC M c {{d u + L B d B )M A - M B d B L A ). (34) 

Equation (|3^ ) introduces numerical problems in the calculation of M A at points where 
J = 0. We avoid these by setting M A = e t,y F A , where 



Substituting F A for M A and setting 5 = a — 7, equation fl3"4[) takes the regularized form 

2i(d u + L A d A )5 = D A L A + H AC F c ((d u + L B d B )F A - F B d B L A ). (36) 
We can now express the "inertial" news ( p0|) in the g a p frame as 

N = h- 2ia u- 2 M a M^( ^ a ^ - uV a Vp^). (37) 

with M a = (0, 0, M A ). An explicit calculation then leads to 

N = ^e~ 2lS uj- 2 e- 2H F A F B {(d u + £ L )c AB + \c AB D c L c + 2ujD A [^ 2 D B {uje 2H )}}, (38) 

where £ l denotes the Lie derivative with respect to L A . The eth versions of these expressions 
are given in Appendix |B[ 
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B. Inertial coordinates 



Equation fl38|) is an expression for the news N(u, x A ) in terms of the angular coordinates 
x A determined by the gauge conditions on the inner worldtube T. The angular coordinates 
y A of inertial observers at X + are constant along the generators of X + , so they are related 
to the x A coordinates by h a d a y A = or 

(d u + L B d B )y A = 0. (39) 

We solve this equation for y A (u,x B ), with the initial condition y A (uo,x B ) = x A , so that y A 
are stereographic angular coordinates. This yields the waveform N'{u, y A ) = N(u, x B (u, y A )) 
measured by inertial observers in terms of the retarded time coordinate u in which the 
evolution is carried out. 

In order to complete the transformation to inertial coordinates we reexpress u in terms 
of Bondi time ub, which is an affine parameter for the generators of X + in a Bondi frame 
satisfying 

fi a d a u B = 1, (40) 

This allows us to rewrite the news as a complex function N(us,y A ) = N'(u(uB,y A ),y A ) of 
inertial (Bondi) coordinates on X + . 

The news can now be readily decomposed into the two standard polarization modes. 
First, we must reexpress the news in terms of a polarization dyad based upon standard 
(6, (ft) spherical coordinates rather than upon stereographic coordinates. The spin rotation 
relating the inertial dyads Q A and (respectively belonging to the south and north stereo- 
graphic patches of the stereographic coordinates y A ) and the standard dyad Q A in spherical 
coordinates (Q a 8a = dg + (i/ sin6*)^) is 

Q A = e^Qj, (41) 
Q A = e^Q A . (42) 

Since the news function has spin- weight 2, the transformation property of spin weighted 
fields [24 1 leads in each patch to the news function M referred to the the dyad Q , 



AT = e - 2i ^N s , (43) 
M = e 2i *N N . (44) 

The decomposition of Af into real and imaginary parts yields the standard polarization 
modes 

N + = U[N] (45) 
iV x =3[A/]. (46) 



V. NUMERICAL IMPLEMENTATION 

In this section we describe a numerical implementation based upon a second-order accu- 
rate finite difference approximation to the equations presented in Sec [Q|. 
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A. The compactified grid 



We introduce a compactified radial coordinate x = r/(R + r), where the scale factor 
R is matched to the radius of the world tube. (In the code runs described here, R is of 
order 1). The radial coordinate is discretized as Xj = X + (i — l)Ax for i = 1 . . . N x , with 
Ax = (1 — X)/(N X — 1). Here setting X = \ corresponds to setting r = R at the inner 
boundary of the grid. The point x^ x = 1 lies at null infinity. The stereographic coordinate 
£ = q + ip is discretized by qj = — 1 + (j — 3) A and p^ = — 1 + [k — 3) A for j, k = 1 . . . N% and 
A = 2/ (iVf — 5). The evolution proceeds with time step Am, subject to a CFL condition 



The fields J, /?, Q and W are represented by their values on this rectangular grid, e.g. 
Jijk = J( u m Xi, QjiPk)- However, for stability [18[] the field U is represented by values at the 
midpoints x i+ i = X{ + Ax/2 on a radially staggered grid (so that U^ k = U(u n , x i+ i, qj,Pk))- 
In the following discussion, it is useful to note that asymptotic flatness implies that the fields 
P(x), U(x), W(x) = r~~ 2 W(x) and J(x) are smooth at future null infinity T + where x — 1. 



B. Hypersurface equations 

We discretize the hypersurface equations (0) - (|ToD along the lines described in |18| 



The only difference here, is the introduction of the nonlinear terms, which are evaluated by 
second order finite differences at the midpoint of the integration cell. 

C. Evolution equation 



The evolution equation (O) has the schematic form 



2 (rJ) >t , - (r- V (r J)_ r ) ^ = ip, + H, (47) 

where we have split off the term 

H = - r -\ r 2 W) r + 2r" 1 e /3 g 2 e /3 - (r^W) r J + Nj- -P x (48) 

r 

which contains no w-derivatives. In order to obtain an accurate global discretization, it is 
convenient to introduce the evolution variable $ = xJ. ($ vanishes in the limit r — > and 
is smooth at X + ). With this substitution, the evolution equation becomes 

2 ($ s (l - x) + $),„ - J^i($ x (l - x) + $) - F 2 q> tXX 
= J(l - x) (-| ^, U (J X K - JK, X )} - 8 ^(1 - x + zW)) + W, (49) 



where 



^ = W" jie x(l - x) + W (50) 
F 2 = (l-x) 2 ((l-x) +xW). (51) 
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The elementary computational cell consists of the grid points (n,i,k,l), (n,i ± l,k,l) 
and (n, i — 2, k, I) on the "old" hypersurface and the points (n + 1, i, k,l), (n + 1, i — 1, k, I) 
and (n + l,i — 2,k,l) on the "new" hypersurface. We center the evaluation at the point 
in + 1/2, i — 1/2, k, I) and approximate the derivatives of $ (to second order). For example, 
we set 



i 



((* ? +1 - ^i 1 ) - - sju + $r+i - (52) 



and 



The evolution to points (n + 1, i, k, I) proceeds in an outward march along the character- 
istics. However, the variable $ appears in the right hand side of Eq. ( fj9|) so that it is not 
convenient to adopt the null parallelogram algorithm used in the quasispherical case |18| 



Instead, we use a Crank-Nicholson [^] scheme to solve the finite difference version of Eq. 
([49|) for . Rather than implementing an implicit solver to obtain we use an iter- 

ative scheme that ensures second order accurate evolution. In addition, the use of 4 points 
on level n in the finite differencing of Q tXU in Eq. ( |52"D introduces some numerical dissipation 
that stabilizes the code even in the regime of very high amplitude fields. 



D. Coordinate transformations 

In order to compute the waveforms N + and iV x measured by inertial observers, we need 
to carry out the transformation from the (u,x A ) frame to the (uB,y A ) frame. We achieve 
this by first transforming to the (u,y A ) and then to (ub,V A ), as follows. 

• Angular transformation: Rather than solving the partial differential equation fl39|) it 
is more expedient to solve the ordinary differential equation 

^ = L>,**) (53) 

with the initial condition x a {uq) = y A and solutions written as x A (u, y A ). This tracks 
a given curve (the characteristic of the PDE) along which y A = const without having 
to track the entire congruence. Since 

dy A = d u y A du + d x sy A dx B 

= du(d u + L B d xB )y A J (54) 
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fl53|) is entirely equivalent to (p^). 

The numerical integration of (|53|) is straightforward and is implemented by the centered 
second order accurate scheme: 

(x* 1 / 1 - x n A ) l = Au (L A )f**\ (55) 



Since we are covering the sphere with two stereographic patches pi |, care must be 
taken to track the y A = const curves as they pass from one patch to the other. This 
is accomplished by the use of a mask which assigns the patch in which the curve lies 
at any time step. 

Bondi time transformation: Similarly, rather than solving the partial differential equa- 
tion (^0[), we determine the transformation to Bondi time by means of the equivalent 
ordinary differential equation du/dus = n a d a u, where the differentiation is along the 
generators of X + . Expressing the left hand side in the (u,£, x A ) coordinate system, 
this reduces to the differential equation du/duB = oJ~ x e~ 2H . We set ub\ u =u — u o as 
the initial condition. 

Again, the integration of this equation is straightforward. We use the centered second 
order accurate algorithm 



-uD^Auiue 2 *)? 1 ^. (56) 



Inertial news: The solution of the above equations determine x A (u, y A ) and u(ub, y A ) 
and thus allow an explicit transformation of the news from grid coordinates (w, x A ) to 
inertial coordinates (uB,y A )- 



VI. NUMERICAL TESTS 
A. Stability 

Numerical experiments were carried out which confirm the stability of the code, subject to 



a CFL condition [18|, in the regime where caustics and horizons do not form. The tests were 



based upon long term evolution of highly nonlinear initial null data with interior worldtube 



data corresponding to the past horizon of a Schwarzschild black hole. See Section |VII| for 
a description of the problem. Stability persists in the regime that the u-direction, along 
which the grid flows, becomes spacelike. This occurs routinely in the evolution of X + , 
on which every tangent direction is spacelike except the null direction of the generators, 
d u + L b 8b- Furthermore, a test evolution of the Schwarzschild spacetime in rigidly rotating 
coordinates remained stable outside the velocity of light cylinder where the grid flow becomes 
superluminal. 
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B. Accuracy 



We have tested the code for accuracy and have verified second order convergence in grid 
size to analytically known values in the following test beds which have been established for 
null codes. 

In all cases, the L2 and norms of the difference between the analytical (F a ) and 
numerical (F n ) solutions was calculated for different grid resolutions. These norms were 
then plotted against Ax and their convergence rate calculated. For every case, a rate of 
at least 1.9 was measured. This procedure is illustrated in the nonlinear axisymmetric 
waveforms case. 

We also checked that the relative error \(F n — F a )/F a \ (when F a 7^ 0) behaved as n(Ax) 2 
with k of order 1, thus establishing the accuracy of the obtained evolution. In particular, 
for the case of Schwarzschild in rotating coordinates, k < 0.7 

• Linearized Waves: 

We used solutions of the linearized characteristic equations on a Minkowski back- 
ground jn|. To check accuracy of the code, we chose a solution representing an outgo- 
ing wave with angular momentum / = 6 with a very small amplitude (| J\ ~ 10 -9 ). The 
inner boundary was set at R — 1. The solution was evolved numerically for Au = 0.5 
and checked to be second order convergent to the analytic linearized solution. 

• Boost and rotation symmetric solutions: 

A family of solutions with exact boost and rotation symmetry called SIMPLE (17 



was also used to check accuracy. Because of their cylindrical symmetry they are 
not asymptotically flat but were used to construct an asymptotically flat solution by 
smoothly pasting to asymptotically flat null data outside some radius R Q . The resulting 
solutions provided a nonlinear test of the second order accuracy of the evolved field 
variables in the domain of dependence interior to R Q . 

Schwarzschild in rotating coordinates: 

By the transformation <ft — > <p + uju of the azimuthal coordinate, the Schwarzschild line 
element in null coordinates becomes 

ds 2 = -(1 uj 2 r 2 sin 2 0)du 2 - Idudr + 2cur 2 sin 2 ddudcf) + r 2 q AB dx A dx B . (57) 

r 



This gauge transformation gives a nontrivial value for U at X + and thus provides a 
useful test bed to check that the numerically calculated wave forms remain zero. We 
evolved this spacetime from u = to u = 0.5 with an inner boundary at R = 3m 
(choosing m = 1 and u — 1). We confirmed that the news calculated at u = 0.5 
converged to zero to second order. This also checked the accuracy of the transformation 
from code coordinates to inertial coordinates at X + . 

• Waveforms from perturbed Schwarzschild black holes: 



14 



The Robinson- Trautman spacetimes describe a distorted black hole emitting purely 
outgoing radiation. The analytic waveform can be easily calculated in the perturba- 



tive regime [IS] and then compared with the results from the code. We evolved the 
spacetime from u = to u = 1, with an inner boundary at R = 3m (where the final 
black hole has mass m — 1). Second order convergence of the numerically obtained 
news function was confirmed under the L\ norm. 

Nonlinear axisymmetric waveforms 

The N x polarization mode must vanish in spacetimes with a twist-free axisymmetry 
(the problem originally studied by Bondi). Since our code uses stereographic coordi- 
nates, axisymmetry is not built in and numerical noise will be produced in the N x 
mode. We then studied the nonlinear scattering of an (£ = 2, m = 0) wave off a 
Schwarzschild black hole and determined to what order this noise converged to zero. 
We gave data at u = and evolved to u = 0.4. The computation was performed on 
grids of size N x — 8n + 1 by N% — 8n + 5 (with n > 5). Figure [I] shows the plot of 
log(|]iV x ||i, a ) vs. log(Ax). Convergence to second order was verified by the slope of 
1.99 

In the next section, we further investigate this scattering problem in the nonaxisym- 
metric case. 



C. Timing 

The characteristic algorithm uses only two complex and two real unknowns (J and U 
and (3 and W, respectively). In addition, the marching structure of the algorithm allows 
us to carry all temporary arrays in two dimensional form. As a result, it is much more 
computationally efficient than Cauchy algorithms for general relativity. Of course, being a 
3-dimensional algorithm, it cannot be run in any practical sense on the current generation of 
workstations or personal computers. On one processor of the C90 machine at the Pittsburgh 
Supercomputing Center, one time step takes 6 seconds for a spatial grid of 64 3 points. A 
physically complete run, such as the ones discussed in the next Section, takes 1 hour (running 
on 1 processor). 



VII. NONLINEAR SCATTERING OFF A SCHWARZSCHILD BLACK HOLE 

The characteristic initial value problem on an outgoing null hypersurface requires an 
inner boundary condition on a worldtube (or, in a limiting case, on the worldline traced out 
by the vertex of a null cone). Some of the worldtube boundary data correspond to gauge 
freedoms. As discussed in Section [TII[ when the worldtube is timelike this is determined by 
the lapse and shift arising from a 2 + 1 decomposition. Here we consider an example in which 
the inner boundary T consists of a ingoing nullcone. We adopt the analogue of setting the 
shift to zero by choosing coordinates x A which follow the ingoing null geodesies. We foliate 
T by slices separated by constant affine parameter A as the analogue of fixing the lapse. We 
then complete the specification of the inner boundary by choosing geometric data on V to 
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correspond to the ingoing r = 2m surface (the past horizon) in a Schwarzschild spacetime. 
This determines inner boundary values satisfying the conservation conditions. 

This construction also determines outgoing null coordinates: u which labels the outgoing 
null hypersurfaces emanating from the the foliation of Y with u\t = A, x A which assigns 
stereographic coordinates to the outgoing null rays and the surface area coordinate r. 

In these coordinates, the Schwarzschild line element takes the Eddington-Finklestein 
form 

(2mA 
1 - — j du 2 - 2dudr + r 2 q AB dx A dx B , (58) 

The initial data corresponds to setting J = as null data on u = 0, with the boundary 
conditions (3 = U A = and W = —2m on T (given by r = 2m). 

We pose the nonlinear problem of gravitational wave scattering off a Schwarzschild black 
hole by retaining these boundary conditions on V (thus continuing to satisfy the conservation 
conditions) but choosing null data at u = corresponding to an incoming pulse with compact 
support, 



All - J 1-— d -— 2 Y l>m r e [R a , R b ] 

j {u = ^ x a )= ) \ rj V rj V 2Z + 1 (59) 



RaV A R b V I 4vr 



otherwise 



where 2^/,m is the spin two spherical harmonic. 

We have calculated the plus and cross components of the news function from u = to 
u = 0.2 (fixing I = m = 2) for various values of A. As expected, for small A, the news 
function scales linearly. However, for larger values, the behavior of the news function shows 
a stronger than linear dependence on A. This is illustrated in Figure |2], which shows the 
plus component of the news, rescaled as N + /X, versus Bondi time. For A < 10~ 2 , the 
rescaled news are essentially equal. In this perturbative regime, the news results from the 
backscattering of the incoming pulse off the effective potential of the interior Schwarzschild 
black hole ||27|| . However, for higher values of A the waveform behaves quite differently. Not 
only is its amplitude greater (up to 3 times for A = 15) but it also reveals the nonlinear 
generation of additional modes at early times for the high amplitude cases. In this regime, 
the mass of the system is dominated by the the incoming pulse, which essentially backscatters 
off itself in a nonlinear way. 

The nonlinearity of the news function can be also seen in Figure |3|, which displays iV x / A 
for amplitudes A = 1.5 x 10~ 6 and A = 15. If the scaling were linear, |iV x /A| would obtain 
a maximum value of ~ 4, while in fact its maximum is ~ 13. 

It is worth emphasizing that the plus and cross components of the news function have 
been obtained by evolving the system in the (w, x A ) coordinate system and then transforming 
to the (ub,V A ) coordinates and to a standard Bondi conformal frame. The surface plot 
of the plus component of the news in the southern hemisphere at Ub = 0.01, given in 
Figure ^, illustrate the smoothness with which our algorithm accomplishes this complete 
transformation. 

ACKNOWLEDGMENTS 



16 



We are grateful to the referee for comments that have improved the clarity of presenta- 
tion. This work has been supported by the Binary Black Hole Grand Challenge Alliance, 
NSF PHY/ASC 9318152 (ARPA supplemented) and by NSF PHY 9510895 and NSF INT 
9515257 to the University of Pittsburgh. N.T.B. and M.M. thank the Foundation for Re- 
search Development, South Africa, for financial support, and the University of Pittsburgh 
for hospitality. J.W. thanks the Universities of South Africa and of Durban- Westville for 
their hospitality. Computer time for this project has been provided by the Pittsburgh Su- 
percomputing Center under grant PHY860023P to J.W. and by a grant from the Center for 
High Performance Computing of the University of Texas at Austin, to Richard Matzner. 

APPENDIX A: NONLINEAR TERMS 

We have the following spin-weighted expressions for the nonlinear terms in the Einstein 
equations: 

R rr = leads to Eq. (|7|) for /? r with 

Nf) = I {JrJr ~ Kl) . (Al) 
o 

R r AQ A = leads to Eq's. (§) and (§) for U ;T and Q, r with 

e 2/3 

Nu = — (KQ-Q- JQ) , (A2) 



and 



N Q = r 2 ^(1 - K)(dK r + 9J r ) + 3(JJ r ) + H{JK r ) - J r M 

+ 2^@J{Jr ~ J 2 lr) + 9J( Jr ~ ^ J,r))J ■ (A3) 

R AB h AB = leads to Eq. (0) for W >r with 

N w = e 2/3 ((1-K) (38/3 + 3,93/3) + X - ( J(8/3) 2 + J(3/3) 2 ) 

~(p/3(M - 3 J) + 3/3(8if - 3 J)) + ^ (j8 2 /3 + .JQ 2 f3 



- e - 2 ^(2KU r U ir + JU 2 + JU 2 r ). (A4) 



RABq A q B leads to Eq. ([12]) for J ur with 



Nj = N JX + N J2 + N J3 + N M + N J5 + N M + N J7 + - (P 1 + P 2 + P 3 + Pa) ( A5) 

r 



where 
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e M /- - -- _ _ \ 

N n = — \K(dJQ0 + 2dKQp - 3J3/3) + J(3 J3/3 - 2QKdp) - JSJS^ 

N J2 = -■= (pj(rU r + 2U) + dJ(rU r + 227)) , 
N J3 = (1 - K)(rdU r + 2W) - J(rm r + 2W), 

N M = r ^e- 2 ?(K 2 Ul + 2JKU r U r + J 2 U 2 r ), 



Nj 5 = - T -J r (dU + W), 

N J6 = r(~{UQJ + UHj)(Jj r - J J, 



+(JK >r - KJ, r )UdJ - £7(9J )P - 2KBKJ r + 2JdKK, r ) 
-U(8j r - K3JJ r + J3 JK yr )j , 

N J7 = r(J tT K - JK^ r ) (7/(8 J - dK) + U(M - 5 J) 
+K(W - W) + ( JW - J9£/)) , 

Pi = r 2 {^{J, r K - JK >r ) + -^{J, r K - JK^ - 8VP >r , 

P 2 = e 2/3 ( - 2K(mp + 8/35/3) - (B/3QK + d(3dK) 

+ (j(8 2 /3 + (8/3) 2 ) + J(3 2 /3 + (5/3) 2 )) + (8 J8/3 + 8 J8/^ , 
p 3 = - ((r8rj r + 2§C/) + (rW r + 29t7)) , 



4 

r -e" 2/3 (2KU r U r + JC/ 2 + JE/J ) . (A6) 
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APPENDIX B: CALCULATING THE NEWS IN THE ETH FORMALISM 

The eth versions of the expressions needed to calculate the news via Eq. (^) are given 
below: 

H AB D A D B log uo = - (-2d 2 log uoJ - 28 2 log uoJ + 488 log uK 

- 8 logo;8 J J 2 — 3 logcuS J J J — 281ogu;3J 

+ 28 log u<5KJK + 8 log ujBJJK + 8 log ujQ JJK 

- 28 log ujBKJJ + 8 J8 log oo J K + 3 J8 log uJK 

- 23/^8 log ujJJ - 8 log ujSJJJ - 28 log cu8 J 

- 81ogcu8JJ 2 + 2d log loBK J K), (Bl) 

D A U A = [W + W)/2, (B2) 
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H AC F c d u F A = (-J U JK - J J + K,JJ + KJC + K :U )/(JJ + 2K + 2). 
For any quantity $ having spin-weight zero, 

u A d A <& = -(um + um). 

The news N is: 

N = ^e~ 2i5 uj- 1 A- 1 {s 1 + s 2 + -(W + W)s 3 - Auj~ 2 s a + 2uj- 1 s 5 ), 

where A = ue 2lS and the terms s 1 to s 5 are: 

si = (J 2 J in + JJJeu - 2JKK /U - 2JK /U + 2J /u K + 2J eu )/(K + 1), 

s 2 = (QJjJJU + 2dJ / KU + 2dJ e U + dJ e J 2 U 

- 2dK : iJKU - 2dK ie JU + 2dUJJK : e - 2dUJJ e K 

- 2WJJ e + AWKK^ + AWK ti + 2dUJJJ / - 2dUJKK / 

- 2WJK / + 4WJ e K + AdUJe + dJ, e JJU + 2dJ e KU + 2§ J e U 
+ dJ e J 2 U - 2dK i£ JKU - 2dK : iJU + 2WJ 2 J^ - 2WJKK yl 

- 2WJK y£ + 2WJ 2 K / - 2WJJ e K - 2dUJJ e )/(2(K + 1)), 

s 3 = (J 2 J,i + JJJ,e ~ 2JKK 4 - 2JK i£ + 2J, e K + 2J e )/(K + 1), 

s 4 = (dAdujJJ + 2dAdujK + 29^9^ - dAdcoJK 

- QAQluJ - dudAJK - QuQAJ + QAQuJ 2 )/(2(K + 1)), 

s 5 = (25 2 A J J + 45 2 Aftf + 4d 2 A + 25 2 A J 2 - 4§5A JK 

- 4SSAJ + J J J 2 + 25A5 JJK + 2<M<3JJ + 9A9JJ 2 J 

+ 29,45 JJfsT + 25A5 JJ - 23A5fT J JK - J J - AdAdKK 

- 4<5A<5K - 3A3JJJK + 2QAQJK + 2<M3J - 3AdJJ 2 K 
+ 2dAdKJ 2 J - dJdAJJK - 2dJdAJJ - 2dJdAK 

- 23JM - dMAJ 2 K - 23JMJ 2 + 2mMJ 2 J 
+ + 4dKMj + MdJJ 2 J 

+ QAQJJ 3 - 2QAQKJ 2 K)/(A(K + 1)). 
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FIGURES 
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FIG. 1. The log of the L2 norm of iV x vs. log(Ax). The computation was performed on grids 
of size N x = 8n + 1 by N% = 8n + 5 (with 5 < n < 13). The corresponding slope is 1.99, indicating 
second order convergence. 
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FIG. 2. Rescaled (+) component of the news for A = 1.5 x 10 n , with n = —2,-1,0,1, and 
A = 7.5 in a 2^22 mode at 6 = (ft = 45°. For the smallest values of A, N + scales linearly with A, 
while at higher values it shows significant nonlinearity which generates additional modes at early 
times. 
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FIG. 3. Rescaled (x) component of the news, for A = 1.5 x 10™, with n = —6 and n = 1 in a 
2Y22 mode at = 4> = 45°. N x /X for A = 15 differs by about 300% from the rescaled value for 
A = 1.5 x 10~ 6 . 
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FIG. 4. Surface plot showing N + (ub = 0.01, y A ) on the southern hemisphere. The obtained 
news function remains smooth after performing the sequence of transformations described in 
Sec. IV B which takes us from a (u, x A ) frame to a Bondi frame (ub,V A )- 
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